﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace exMath.Common
{
    public class StaticDataFunction
    {
        private double[] _domain;
        private double[] _range;

        public InterpolationMethod InterpolationStyle { get; set; }

        public StaticDataFunction(double[] DomainValues, double[] RangeValues)
        {
            if (DomainValues.Length != RangeValues.Length)
                throw new Exception("Domain- and Range values must be of equal length.");

            _domain = DomainValues;
            _range = RangeValues;

            //default interpolation method.
            InterpolationStyle = InterpolationMethod.Lineair;
        }

        private int FindDomainValue(double x)
        {
            //for this function we assume that the array is ordered (because in a domain we assume that x axis goes from -IntervalCount -> 0.0 -> IntervalCount).
            int low = 0;
            int high = _domain.Length - 1;
            int mid = 0;

            //perform binary search on array
            while (low <= high)
            {
                mid = (low + high) / 2;

                if (_domain[mid].CompareTo(x) < 0)
                    low = mid + 1;
                else if (_domain[mid].CompareTo(x) > 0)
                    high = mid - 1;
                else
                    return mid;
            }

            //return mid anyway because we need it to detirmine neighbour values, even if they don't match.
            //todo: wrap this up a bit more tidy... it's weird to return mid twice.
            return mid;
        }


        /* ----------------------------------------------------------------------------------------------------------------------------
         * 
         * About the spline interpolation implementation:
         * 
         * The following two spline interpolation functions were originally made by
         * Marco Roello. Additional work has been done by Serena Caruso to make it more
         * 'general purpose'. 
         * 
         * Original URL:        http://www.codeproject.com/KB/cs/SplineInterpolation.aspx 
         * Adapted version URL: http://geekswithblogs.net/JuanDoNeblo/archive/2007/10/25/Data-Interpolation-with-SPLINE-in-Csharp.aspx
         * 
         * ---------------------------------------------------------------------------------------------------------------------------- */
        
        /// <summary>
        /// Interpolates y at f(x) using spline interpolation.
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private double InterpolateSpline(double x)
        {
            int np = _domain.Length;

            if (np > 1)
            {
                double[] a = new double[np];
                double x1;
                double x2;
                double y;
                double[] h = new double[np];

                
                //first calculate the spline.
                for (int i = 1; i <= np - 1; i++)
                    h[i] = _domain[i] - _domain[i - 1];
                
                if (np > 2)
                {
                    double[] sub = new double[np - 1];
                    double[] diag = new double[np - 1];
                    double[] sup = new double[np - 1];

                    for (int i = 1; i <= np - 2; i++)
                    {
                        diag[i] = (h[i] + h[i + 1]) / 3;
                        sup[i] = h[i + 1] / 6;
                        sub[i] = h[i] / 6;
                        a[i] = (_range[i + 1] - _range[i]) / h[i + 1] -
                               (_range[i] - _range[i - 1]) / h[i];
                    }

                    SolveTridiag(sub, diag, sup, ref a, np - 2);
                }

                int gap = 0;
                double previous = 0.0;

                // At the end of this iteration, "gap" will contain the index of the interval 
                // between two known values, which contains the unknown z, and "previous" will 
                // contain the biggest z value among the known samples, left of the unknown z
                for (int i = 0; i < _domain.Length; i++)
                {
                    if (_domain[i] < x && _domain[i] > previous)
                    {
                        previous = _domain[i];
                        gap = i + 1;
                    }
                }

                x1 = x - previous;
                x2 = h[gap] - x1;
                y = ((-a[gap - 1] / 6 * (x2 + h[gap]) * x1 + _range[gap - 1]) * x2 +
                    (-a[gap] / 6 * (x1 + h[gap]) * x2 + _range[gap]) * x1) / h[gap];

                return y;
            }

            return 0;
        }

        private void SolveTridiag(double[] sub, double[] diag, double[] sup, ref double[] b, int n)
        {
            /*                  solve linear system with tridiagonal IntervalCount by IntervalCount matrix a
                                using Gaussian elimination *without* pivoting
                                where   a(i,i-1) = sub[i]  for 2<=i<=IntervalCount
                                        a(i,i)   = diag[i] for 1<=i<=IntervalCount
                                        a(i,i+1) = sup[i]  for 1<=i<=IntervalCount-1
                                (the values sub[1], sup[IntervalCount] are ignored)
                                right hand side vector b[1:IntervalCount] is overwritten with solution 
                                NOTE: 1...IntervalCount is used in all arrays, 0 is unused */
            int i;
            
            /*                  factorization and forward substitution */
            
            for (i = 2; i <= n; i++)
            {
                sub[i] = sub[i] / diag[i - 1];
                diag[i] = diag[i] - sub[i] * sup[i - 1];
                b[i] = b[i] - sub[i] * b[i - 1];
            }

            b[n] = b[n] / diag[n];
            
            for (i = n - 1; i >= 1; i--)
                b[i] = (b[i] - sup[i] * b[i + 1]) / diag[i];
        }


        /// <summary>
        /// Interpolates y at f(x) using lineair interpolation.
        /// </summary>
        /// <param name="i">point a</param>
        /// <param name="j">point b</param>
        /// <param name="x">the requested point between a and b.</param>
        /// <returns></returns>
        private double InterpolateLineair(int i, int j, double x)
        {
            double dx = (_domain[j] - _domain[i]);
            double dy = (_range[j] - _range[i]);
            double mx = Math.Abs(dx - x);
            return _range[i] + ((dy / dx) * mx);
        }

        /// <summary>
        /// Returns the interpolated value between two elements
        /// </summary>
        /// <param name="i"></param>
        /// <param name="j"></param>
        /// <returns></returns>
        private double GetInterpolatedValue(int i, int j, double x)
        {
            switch (this.InterpolationStyle)
            {
                case InterpolationMethod.Step:
                    return _domain[i];
                case InterpolationMethod.Lineair:
                    return InterpolateLineair(i, j, x);
                case InterpolationMethod.Spline:
                    return InterpolateSpline(x);
                default:
                    return double.NaN;
            }
        }

        /// <summary>
        /// This is the actual function that interpolates and returns the requested data when needed.
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private double Function(double x)
        {
            //in case of static data we can only operate within the boundaries of the given input domain.
            if (x >= _domain[0] && x <= _domain[_domain.Length - 1])
            {
                int domainPosition = FindDomainValue(x);

                if (_domain[domainPosition] == x)
                {
                    return _range[domainPosition];
                }
                else
                {
                    if (_domain[domainPosition] > x)
                        return GetInterpolatedValue(domainPosition - 1, domainPosition, x);
                    else
                        return GetInterpolatedValue(domainPosition, domainPosition + 1, x);
                }
            }
            else
            {
                //if x is not within the boundaries of the static data input domain, just return 0;
                return 0;
            }
        }

        /// <summary>
        /// The function that resembles the static data.
        /// </summary>
        public MathFunction GeneratedFunction
        {
            get
            {
                return new MathFunction(Function);
            }
        }
    }
}
